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Multi-chamber ,  lumped  parameter  models  of  the  pulmonary  and 
circulatory  systems  have  been  developed  which  incorporate  effects  of 
whole-body  accelerations  (WBA)  for  any  body  orientation  and  a  variety 
of  breathing  maneuvers.  The  pulmonary  model  has  been  used  in  an 
extensive  study  of  effects  of  sustained  WBA  on  pulmonary  mechanics  and 
ventilation.  This  study  Investigated  the  effectiveness  of  altering 
seat  back  angle  to  Increase  G-tolerance. 

The  development  snd  testing  of  the  separate  pulmonary  and  circu¬ 
latory  models  has  been  reported  in  previous^AFOSR* reports,  as  well  as 
the  above  described  WBA  study. 

This  final  report  describes  latest  efforts  to  incorporate  feed¬ 
back  and  control  mechanisms  in  the  circulatory  model  to  simulate  the 
baroreceptor  and  chemoreceptor  actions  regulating  blood  flow  and  blood 
pressure  in  the  human  body.  This  work  is  not  yet  complete. 

Also  described  is  a  partially  completed  study  using  the  pulmonary 
model  to  test  and  implement  a  non- invasive  technique  for  diagnosing 
ventilation  distribution  dysfunctions.  ^ 
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PART  I 

FEEDBACK  AND  CONTROL  IN  THE  CARDIOPULMONARY  SIMULATOR 
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INTRODUCTION 

A  computer  simulator  of  the  human  cardlo-pulmonary  system  has 
been  developed  as  two  distinct  models.  The  first  model  Is  a  pulmonary 
simulator  which  has  proved  to  be  a  valuable  stand-alone  system  for  the 
evaluation  of  pulmonary  mechanics  and  ventilation  during  whole  body 
acceleration  (NBA).  This  pulmonary  simulator  has  been  described  In 
detail  in  previous  reports  and  a  working  copy  has  been  sent  to 
Wrlght-Patterson  AFB.  The  second  model  ,1s  a  circulatory  simulator 
which  has  also  been  described  In  detail  In  previous  reports. 

The  achievements  during  this  report  period  are  1)  modification  of 
the  circulatory  model  to  Increase  Its  controllability.  2)  the  Inclu¬ 
sion  of  feedback  mechanisms  to  the  circulatory  system,  and  3)  in¬ 
terfacing  the  pulmonary  and  circulatory  models  to  work  in  conjunction 
with  each  other  as  a  single  simulator. 

A  comprehensive  numerical  model  of  the  feedback  and  control 
systems  which  regulate  the  cardiovascular  system  was  devised  and 
incorporated  into  the  digital  computer  simulator  of  the  cardiovascular 
and  pulmonary  systems.  The  model  is  designed  specifically  to  study 
the  effects  of  short  term  whole-body  acceleration  on  the  cardiovascu¬ 
lar  system  as  encountered  by  astronauts  and  pilots  of  high  performance 
aircraft.  The  feedback  systems  include  the  sympathetic  and  parasympa¬ 
thetic  neural  control  of  the  cardiovascular  system,  resulting  from 
feedback  loop  responses.  These  feedback  loop  responses  are  due  to 
1)  arterial  pressure  changes  acting  upon  specialised  sensors,  i.e. 
baroreceptors ,  and  2)  changes  in  blood  gas  concentrations  acting  upon 
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cheaoreceptors.  The  controlling  mechanisms  include  variations  of 
heart  rate,  contractility,  and  the  peripheral  resistances  and  compli¬ 
ances  in  the  vascular  system. 

The  completed  cardio-pulmonary  computer  simulator  is  a  multi- 
chambered  model  representing  the  human  cardiovascular  and  pulmonary 
systems,  including  pulsatile  blood  flow  and  gas  transport  and  ex¬ 
change.  It  represents  an  active,  time  varying,  nonlinear  model  with 
non-linear  elastic  characteristics,  pressure  dependent  flow  resis¬ 
tances,  and  arterial  pressure  and  blood  gas  concentration  feedback 
systems. 

BACKGROUND 

Many  models  have  been  developed  to  represent  the  pulmonary  or 
cardiovascular  systems;  however,  none  address  the  same  functional 
parameters  or  level  of  complexity  as  the  current  model.  Some  models 
are  limited  to  specific  anatomical  features,  such  as  heart,  systemic 
arterial  vasculature,  or  venous  vasculature.  Other  models  are  de¬ 
signed  to  study  different  functions,  such  as  fluid  balance,  long  term 
regulation,  or  metabolism.  Few  models  include  the  complexity  of 
multi-elements ,  nonlinear  relationships,  or  dependent  elements.  Of 
the  models  including  gravitational  effects  on  the  cardiovascular 
system,  none  have  been  as  complete  as  the  current  model. 

Detailed  descriptions  of  the  pulmonary  model  and  the  cardiovascu¬ 
lar  model  have  been  reported  previously.  This  report  will  review  some 
of  these  details  to  facilitate  describing  the  Inclusion  into  the 
circulatory  model  of  feedback  and  controllability  systems. 
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PHYSIOLOGICAL  BACKGROUND 

A  cardiovascular  feedback  system  can  be  modeled  as '  a  controlled 
system  (plant),  a  feedback  loop,  and  a  controlling  mechanism.  The 
controlled  system  represents  that  portion  of  the  cardiovascular  system 
to  be  regulated.  The  feedback  loop  portion  consists  of  a  receptor 
stimulus  (a  parameter  from  the  controlled  system),  feedback  receptors, 
receptor  response,  any  modulating  effects,  and  the  resultant  stimulus 
to  activate  the  control  mechanism.  The  control  mechanisms  are  the 
physiological  parameters  which  are  altered,  due  to  the  control  mecha¬ 
nism  stimuli,  resulting  in  the  regulation  of  the  controlled  system. 
Figure  1  shows  a  block  diagram  of  a  feedback  system. 

There  are  many  cardiovascular  feedback  systems.  However,  in 
limiting  the  model  to  the  study  of  short  term  effects  due  to  whole 
body  acceleration,  two  feedback  systems  exhibit  the  predominant 
effects.  They  are  the  baroreceptor  reflex  and  the  chemoreceptor 
reflex  (see  Figure  2). 

The  baroreceptor  reflex  is  particularly  critical  in  regulating 
gravitational  or  orientational  changes.  This  feedback  system  is 
stimulated  by  arterial  pressure  and  results  in  controlling:  1)  heart 
rate,  2)  cardiac  contractility,  and  3)  vascular  resistance  and  compli¬ 
ance. 

The  arterial  pressure  in  the  carotid  sinus  and  aortic  arch  is  the 
stimulus  acting  on  the  pressoreceptors  (baroreceptors)  located  in  the 
same  areas.  The  efferent  nerve  response  from  the  baroreceptors 
decreases  in  relation  to  a  carotid  pressure  decrease  (Varner,  1958). 
This  receptor  response  then  enters  the  central  nervous  system,  CNS, 
where  it  is  modulated  with  other  CNS  stimuli,  resulting  in  new 
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efferent  neural  signals  to  the  controlling  mechanisms.  For  these 
efferent  neural  signals ,  the  sympathetic  frequency  of  stimulation  to 
the  controlling  mechanism  is  Increased  while  the  parasympathetic 
frequency  of  stimulation  is  decreased.  The  relationship  describing 
the  CHS  Input/output  and  modulation  with  other  efferent  signals  is  not 
clearly  understood.  However,  relations  between  carotid  pressure  and 
sympathetic  firing  (Kezdl  and  Geller,  1968)  and  carotid  pressure  and 
vagal  firing  (Katona  and  Barnett,  1970)  have  been  made. 

The  responses  of  the  three  controlling  mechanisms  due  to  the 
efferent  neural  signals  complete  the  description  of  this  feedback 
system.  The  heart  rate  Increases  in  response  to  Increased  sympathetic 
activity  and  decreased  parasympathetic  activity  (Warner  and  Cox, 
1962).  The  two  neural  paths  have  different  time  delays  before  their 
effects  are  seen.  The  sympathetic  response  can  be  delayed  20  seconds 
while  the  parasympathetic  response  is  immediate  (Wang  and  Borison, 
1947).  The  carotid  pressure  change  has  been  related  to  the  heart  rate 
by  Scher  and  Young  (1970)  and  Glick  and  Braunwald  (1965). 

The  second  controlling  mechanism  is  heart  contractility,  or  force 
of  contraction,  that  is  not  due  to  the  Frank- Star  ling  effect  of 
stretch.  The  contractility  will  increase  with  Increased  sympathetic 
activity  and  decreased  parasympathetic  activity  (Levy,  1963  and 
Rushmer,  1962).  Additionally,  the  contractility  will  vary  due  solely 
to  changes  in  heart  rate  (Koch-Weser  and  Blinks,  1963). 

The  third  controlling  mechanism  of  the  baroreceptor  reflex 


feedback  system  is  vl  sular  r  /lienee,  which  is  proportional  to 
vascular  rasistence.  Thw  vascular  compliance  decreases  resulting  in 
an  increased  vascular  resistance  when  sympathetic  activity  to  vascular 
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tissue  Increases  and  parasympathetic  activity  decreases.  Not  all 
vascular  compliances  change  similarly.  The  coronary  arteries  and 
cerebral  vessels  do  not  constrict  when  the  rest  of  the  vasculature 
does  (Gregg  and  Fischer,  1963,  and  Rapela  and  Green,  1964). 

Another  Important  type  of  feedback  system  is  the  chemoreceptor 
reflex.  It  is  stimulated  by  the  blood  gases  and  results  in  control¬ 
ling:  1)  heart  rate,  2)  contractility,  3)  vascular  resistance,  and 
4)  respiration.  This  blood  gas  feedback  system  can  be  divided  into 
two  loops.  In  the  first,  blood  gases  in  the  aortic  arch  and  carotid 
sinus  are  the  stimuli  of  chemoreceptors  located  in  the  aortic  bodies 
and  carotid  bodies.  A  decrease  in  0^,  an  increase  in  CO^,  or  decrease 
in  pH  will  cause  the  efferent  nerve  activity  from  the  chemoreceptors 
to  the  CNS  to  increase.  As  with  the  baroreceptors ,  the  CNS  modulates 
this  receptor  response  with  other  CNS  stimuli.  With  the  above  gas 
changes,  the  efferent  neural  signals  from  the  CNS  to  the  controlling 
mechanisms  are  altered  such  that  the  sympathetic  frequency  of  stimu¬ 
lation  is  Increased  while  the  parasympathetic  frequency  of  stimulation 
is  decreased  (Von  Euler,  et  al.,  1939). 

The  second  feedback  loop  of  this  blood  gas  feedback  system  has 
its  chemoreceptors  .located  in  the  medulla.  Similar  to  the  loop 
described  above,  an  Increase  in  CO^,  or  decrease  in  0,  or  pH  will 
cause  the  CNS  to  alter  the  frequency  of  efferent  neural  signals  to  the 
controlling  mechanisms.  The  sympathetic  stimuli  Increase  in  frequency 
while  the  parasympathetic  decrease. 

Although  these  two  loops  are  activated  by  the  same  stimuli,  i.e., 
blood  gas  concentration,  they  are  not  equally  as  sensitive  to  the 
individual  gaseB.  Similarly,  even  though  these  two  loops  both  affect 
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the  same  controlling  mechanisms  (heart  rate,  contractility,  vascular 
resistance,  respiration  rate,  and  tidal  volume),  the  degree  to  which 
each  loop  affects  these  control  mechanisms  is  different. 

Basically,  the  effect  on  these  control  mechanisms  is  the  same  as 
described  for  the  baroreceptors.  With  an  increase  in  frequency  of 
sympathetic  activity,  and  a  decrease  in  frequency  of  parasympathetic 
activity,  heart  rate,  contractility,  and  peripheral  resistance  all 
increase.  In  addition,  these  loops  affect  the  respiration.  With  a 
similar  change  in  sympathetic  and  parasympathetic  activity,  the 
respiration  rate  and  the  tidal  volume  increase. 

Many  additional  feedback  systems  exist  in  the  circulatory  system, 
including  autoregulation  and  the  ischemic  reflex.  However,  the  two 
described  above  exhibit  the  predominant  effects  in  the  types  of 
situations  to  be  modeled. 

MODEL  DESCRIPTION 

In  order  to  understand  how  these  feedback  systems  are  modeled,  it 
is  necessary  to  first  describe  the  cardiovascular  model  into  which 
they  are  incorporated.  The  cardiovascular  system  is  represented  by  a 
twenty-four  chamber  model  Including  a  four  chamber  heart  (see  Figure 
3) .  The  chambers  represent  various  portions  of  the  cardiovascular 
system,  including  the  arterial  vascular  network,  several  capillary 
beds,  several  venous  vascular  networks,  and  the  four  chambers  of  the 
heart . 

The  mechanics  of  blood  circulation  for  a  generalized  chamber 
element  can  be  described  using  three  equations.  The  first  treats  the 
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axial  notion  of  the  blood  through  the  blood  vessels.  The  second 
treats  the  transverse  motion  of  the  blood,  vessel  wall,  and  surround¬ 
ing  tissue,  and  the  third  assumes  conservation  of  mass  for  the  blood 
in  the  chambers. 

Axial  Motion 

The  equations  for  axial  motion  of  the  blood  through  a  circulatory 

element  are  based  on  the  conservation  of  momentum  written  in  terms  of 

forces  as  source  terms  (see  Figure  4) .  The  rate  of  change  in  axial 

momentum  of  blood  in  the  vessel,  AlEkSl  can  be  described  in  terms  of 

at 

the  density  of  blood, p  ,  the  vessel  length,  L,  and  the  mean  volumetric 

flow  rate  through  the  vessel,  q.  According  to  Newton's  second  law  of 

motion,  this  change  in  momentum  is ’equal  to  the  algebraic  sum  of  the 

axial  forces  acting  on  the  blood  in  the  vessel.  The  forces  acting 

consist  of  the  pressure  force  of  the  fluid  at  the  inlet,  ,  and  the 

outlet,  the  *orce  gravity,  pgCbj-h^A,  the  inertial  force  of 

1  2 

the  blood  at  the  Inlet,  pq^  ,  and  at  the  outlet,  -pq2  and  the 

viscous  drag  force  at  the  vessel  wslls,  fA  .  Here,  A  is  the  total 

S  8 

surface  area  of  the  vessel;  A  is  the  average  cross  sectional  area  of 

the  vessel;  (h^-h^).  is  the  distance  separating  the  inlet  and  outlet 

ends  of  the  vessel  along  the  direction  of  the  gravitational  accel- 

dx! 

erstlon,  g;  q^  and  q^  ere  inlet  and  outlet  volumetric  flow  rates, 
dx2 

and  are  the  inlet  and  outlet  velocities;  and  f  is  the  total  drag 
force  per  unit  area  of  vessel  wall.  Assuming  a  uniform  cross  section¬ 
al  area,  i.e.,  A^-A^-A,  and  assuming  fA&  can  be  approximated  by  a 
resistance  times  the  mean  volumetric  flow  rate,  q,  the  equation  can 


be  written  as: 


CHAMBER  FLOW  DYNAMICS 


AXIAL  MOTION 
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■  Ml  sr  -  0fl2  p-  *  v  -  v  *  c«(hi '  VA  -  V 


This  can  be  rearranged  Co  get: 


.  .  dx.  2  dx,  2 

A~3?“pdr  ■  p  dT  +  ?!  -  P2  +Pg(h1  -  h2)  -  Rfq 


where  -  R^/A  and  R^  is  a  function  of  volume. 


Transverse  Motion 

The  equations  for  transverse  motion  of  the  vessel  wall  and 

surrounding  tissue  are  also  derived  from  Newton's  second  law  (see 

Figure  5).  The  sum  of  the  forces  acting  on  the  vessel  wall  and 

surrounding  tissue  are  equal  to  the  mass  of  the  wall  and  tissue,  m, 

times  its  acceleration,  a.  In  terms  of  volume  of  the  chamber  element , 

2 

V,  the  term  ma  can  be  expressed  as  —  — •>»  where  A  is  the  surface  area. 

A  dt 

The  forces  acting  on  the  vessel  wall  consist  of  the  force  due  to 
blood  pressure  in  the  vessel,  PA;  the  force  due  to  external  pressure 
on  the  vessel  and  surrounding  tissues,  AP',  the  force  due  to  the 
elastance  pressure  from  the  muscular  tissue  surrounding  the  vessel, 
-K(V) (V-Veq) ;  the  force  due  to  the  friction  of  the  wall  on  surround¬ 
ing  tissues,  -R  4r»  and  in  the  case  of  the  heart  chambers,  an  external 
at 

force  due  to  heart  contraction,  FF(t).  This  gives  the  equation: 


j  —•  -  AP  -  AP'  -  K(V)  (V-Veq)  -  R  ~  -  FF(t) 


(2) 
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Here  K(V)  Is  a  nonlinear  compliance  term  dependent  on  volume;  R  Is  a 
resistance  term  associated  with  frictional  resistance  of  tissue;  and 
FF(t)  is  a  time  varying  function  representing  the  force  due  to  con¬ 
traction  of  the  heart. 

Conservation  of  Mass 

Based  on  conservation  of  mass,  the  rate  of  change  in  volume  of  a 
dV 

chamber,  must  equal  the  difference  between  the  sum  of  the  volu¬ 
metric  flows  into  the  chamber,  q.  ,  and  the  sum  of  the  volumetric 

in 

flows  out  of  the  chamber,  qout«  This  gives  the  equation: 

dV 

-  1  “in  -  E  ’out  (3> 

The  above  equations  represent  a  generalized  chamber  in  the 
twenty-four  chamber  cardiovascular  model.  However,  each  individual 
chamber  can  be  simplified  according  to  the  particular  dynamics  of  that 
chamber.  For  example,  only  the  four  heart  chambers  have  the  forcing 
function  term  FF(t)  representing  an  external  force  due  to  heart 
contraction. 

The  family  of  equations  describing  the  cardiovascular  system  can 
be  illustrated  as  a  simple  block  diagram  with  each  block  representing 
a  particular  portion  of  the  system  (see  Figure  6).  From  such  a  block 
diagram,  it  can  be  clearly  illustrated  where  the  feedback  systems  will 
be  included  (see  Figure  7).  The  dotted  lines  represent  stimuli  to  the 
feedback  receptors.  The  blocks  labeled  FF(t)  represent  the  feedback 
responses  changing  the  forcing  function  to  alter  heart  rate  and 
contractility.  The  blocks  labeled  K(V)  represent  the  feedback 
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MODEL  DESCRIPTION  OF  CONTROL  AND  FEEDBACK  SYSTEMS 

The  physiological  feedback  reflexes  described  earlier  -sense  three 
parameters.  The  carotid  blood  pressure  is  sensed  by  the  baroreceptor 
reflex  and  the  partial  pressure  of  carbon  dioxide,  P  ,  and  the 

partial  pressure  of  oxygen,  P  ,  are  sensed  by  the  chemoreceptors. 

°2 

Both  of  these  reflexes  result  in  controlling  heart  rate,  heart  con¬ 
tractility,  and  vascular  compliance  which  affects  the  vascular  resis¬ 
tance  (see  Figure  8). 

The  model  can  simulate  these  changes  by  varying  just  two  parame¬ 
ters.  first,  the  model  variable  for  the  forcing  function,  FF(t),  of 
each  heart  chamber  can  be  altered  to  control  the  heart  rate  and 
contractility  of  the  heart.  Second,  the  vascular  compliance  term, 
K(V),  of  each  chamber  can  be  altered  to  control  vascular  compliance 
which  will  result  in  a  change  in  vascular  resistance. 

To  model  the  physiological  reflexes,  both  the  forcing  function 
terms  and  the  chamber  vascular  compliance  terms  should  be  appropriate¬ 
ly  varied  in  response  to  changes  in  either  carotid  blood  pressure  or 
the  partial  pressure  of  CO^  or  the  partial  pressure  of  0^  (see  Figure 
7). 

To  simulate  a  four  chamber  pulsating  heart  with  variable  heart 
rate  and  variable  contractility,  a  new  heart  model  was  developed  which 
Included  a  new  forcing  function.  The  forcing  function,  FF(t),  as 
shown  in  Equation  2,  is  a  time  varying  function  representing  the  force 
and  duration  of  a  contraction  of  the  heart.  For  each  different  heart 
rate  and  contractility,  a  unique  forcing  function  is  required. 

The  duration  of  ventricular  contraction  closely  corresponds  to 
the  Q-T  Interval,  Q?lnt*  in  an  electrocardiogram.  An  equation 
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Figure  8 
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relating  heart  rate,  HR,  to  ventricular  contraction  time  was  developed 
using  data  on  Q-T  intervals  (Ashman,  1956). 

Since  the  heart  rate  is  inversely  proportional  to  its  pulse 
interval,  RR^nt>  we  can  write  the  equation  for  the  ventricles  as: 


1.6 


(4) 


where  HR  -  60  x  RR.  „ . 

lnt 

Similar  data  does  not  exist  for  the  duration  of  atrial  con¬ 
traction,  PR^nt ,  but  using  known  values  at  several  heart  rates  and 
assuming  a  slope  similar  to  that  for  the  ventricles,  the  following 
equation  was  derived: 


5.4  ♦ 


RR 


lnt 


(5) 


Now,  for  a  given  heart  rate  the  appropriate  forcing  function,  FF, 
for  each  heart  chamber  can  be  computed. 


FF  -  PSAMP 


•  AlPb.  . 

1  t  -  PSCENT  4U  f 

V.  1  PSEND  -  PSCENT I  J 

«m^  ^ 


PPAVE 


(6) 


where:  t  ■  time 

PSEND  •  time  at  the  end  of  FF,  or  time  at  start  of 

contraction  +  the  duration  of  the  contraction 

PSCENT  -  time  midway  into  FF 

Alpha  -  a  constant  controlling  the  shape  of  the  curve 

PSAMP  ■  variable  controlling  the  amplitude  of  the  FF 
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FPAVE  -  addition  force  on  chamber  due  to  plural  cavity 
pressure. 

Changing  variable  PSAMP  corresponds  to  changing  the  contractility 
of  the  heart  chambers. 

Physiologically,  contractility  or  inotropic  effect  cannot  be 

measured  directly;  however,  a  good  Indication  of  contractility  is  the 

dP 

maximum  rate  of  change  in  the  pressure,  —  Max,  (Guyton,  1971).  By 
infusion  of  isoprenaline,  which  Increases  contractility,  and  by 
controlling  the  heart  rate,  a  linear  relationship  can  be  seen  between 
increase  in  heart  rate  and  4^  Max  (Fumival,  et  al.,  1970).  This 
relationship  was  used  initially  to  modify  PSAMP  to  generate  a  Max 
at  a  given  heart  rate.  However,  it  was  modified  by  making  repeated 
computer  runs  so  that  stroke  volumes  and  cardiac  outputs  would  corre¬ 
spond  to  known  values  under  various  degrees  of  sympathetic  and  para¬ 
sympathetic  stimulation  (Guyton,  1971).  Therefore,  for  a  desired 
heart  rate,  using  equations  4  and  5,  the  duration  of  the  forcing 
functions  and  heart  rate  intervals  can  be  calculated  for  the  four 
heart  chambers,  and  the  degree  of  contractility  for  each  can  be  set  by 
calculating  the  values  of  PSAMP.  Selection  of  the  appropriate  vascu¬ 
lar  compliance  will  be  discussed  later. 

Physiologically,  regardless  of  the  type  of  receptor  stimulus  used 
as  input  to  the  feedback  loop  (blood  pressure,  P__  or  P_  )  the  output 
to  the  controlling  mechanism  will  always  be  in  terms  of  sympathetic 
and  parasympathetic  neural  activity.  Also,  regardless  of  the  type  of 
controlling  mechanism  (heart  rate,  contractility,  or  vascular  compli¬ 
ance)  ,  the  input  stimulus  will  always  be  expressed  in  terms  of  sympa¬ 
thetic  and  parasympathetic  neural  activity. 


This  relationship  implies  that  if  the  resultant  sympathetic  and 
parasympathetic  stimulus  due  to  combined  stimulation  by  carotid  blood 
pressure,  P^,  ,  and  P^  ,  could  be  found,  then  the  controlling  mecha- 
nisms  could  each  be  calculated  from  that.  Much  data  exists  comparing 
sympathetic  and  parasympathetic  activity  to  heart  rate  and  other 
functions.  Also,  there  are  many  investigations  concerning  activation 
of  receptor  sites  and  the  resulting  neural  responses.  However, 
sufficient  data  does  not  yet  exist  to  bring  these  systems,  from 
receptor  to  controlling  mechanism,  into  a  unified  theory  based  on 
sympathetic  and  parasympathetic  activity.  A  primary  problem  is  a  lack 
of  quantitative  information  on  Interactions  in  the  central  nervous 
system. 

The  approach  used  in  this  model  was  to  take  a  lumped  systems 
approach  by  looking  at  the  receptors  stimulus  (carotid  blood  pressure) 
and  a  controlling  mechanism  output  (heart  rate).  Initially,  the  mean 
carotid  pressure,  P£  ,  at  the  baroreceptors  is  calculated  for  the 
current  heartbeat.  The  carotid  pressure  is  calculated  as  the  arterial 
chamber  pressure  plus  the  hydrostatic  pressure  between  the  heart  and 
the  carotid  artery  where  the  baroreceptors  are  located.  This  pressure 
is  compared  to  a  set  normal  functioning  carotid  pressure,  P^  .  The 
difference  in  carotid  pressure,  APc»  is  used  to  compute  the  change  in 
heart  rate  necessary  to  return  the  carotid  pressure  to  its  normal 
value.  The  change  in  heart  rate,  A  HR,  is  the  difference  between  the 
normal  operating  heart  rate,  HRj  ,  and  the  new  heart  rate,  HR,  needed 
to  correct  the  pressure. 


23 


The  following  equations  give  results  similar  to  existing  data  on 
change  in  carotid  pressure  versus  change  in  heart  rate  (Glick,  1965 
and  Robinson,  1966). 


For  AP  >0 
c  — 


A  HR  -  e 


-a(APc) 


-  1  - 


<APC> 


(7) 


For  AP  <0 
c 


AHR  =  e 


-a(APc) 


-  1 


Since  changes  in  the  partial  pressure  of  C02  and  0^  also  result 

in  changes  in  heart  rate,  the  equation  must  be  amended  to  reflect 

this.  For  the  blood  gases  and  chemoreceptor  reflex,  similar  normal 

operating  values  are  set  for  the  partial  pressure  of  C02»  P^  j,  and 

the  partial  pressure  of  0^,  P^  j  .  Subtracting  these  from  their 

respective  current  values  gives  the  change  in  partial  pressure  of  0^, 

AP_  and  the  change  in  partial  pressure  of  C02,  APqq  .  Including  the 
°2  2 
blood  gas  effects  into  the  equation  above  results  in  the  following 

equation: 


1 

> 


-aAP 


AHR 


'  -  1  -  4P=  ♦  CiPC02  -  d4P02 


(8) 


As  described  above,  once  the  desired  heart  rate  is  known,  then 
the  forcing  function  can  be  adjusted  to  provide  the  proper  heart 
durations  and  contractility. 


’AVlV,  v  - 


.j. 
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The  only  remaining  controlling  mechanism  is  the  vascular  compli¬ 
ance,  K(V) .  After  the  other  control  mechanisms  were  changed  accord¬ 
ingly,  the  vascular  compliances  and  their  rate  of  variability  were 
found  by  making  multiple  computer  simulations. 


CONCLUSIONS 

The  pulmonary  simulator  has  proved  to  work  independently  in 
evaluation  of  pulmonary  mechanics  and  ventilation  during  whole  body 
acceleration.  The  circulatory  model  and  pulmonary  model  work  in 
conjunction  with  each  other  as  a  single  simulator.  Inclusion  of 
feedback  mechanisms  into  the  circulatory  model  improves  its  perfor¬ 
mance  . 
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APPENDIX  TO  PART  I 


This  appendix  lists  the  subroutines  and  their  purposes  followed 
by  a  complete  computer  listing  of  the  programs.  The  outlining  of  the 
subroutines  will  not  include  the  details  of  the  lung  program.  A  flow 
chart  describing  the  order  of  calculation  is  shown  in  Figure  27. 
BMAIN  -  program  of  pumping  heart  and  blood  flow  only 
(circulatory  portion  of  the  model  only) 

BMAIN  -  calling  routine  controlling  storage,  subroutines,  time 
incrementation,  and  1.0.  Subroutines  of  BMAIN 

BDRIVS  -  calculates  flows,  derivatives  of  volumes,  and  controls 
the  leg  valve. 

RESIST  -  calculates  resistance  of  all  segments.  Options  avail- 

2 

able  are  LOPR*neg  makes  R  •  k/(V  +V).  This  is  to  model  non- 
Newtonian  flow  in  the  capillaries  if  desired.  Not  presently  in 
use. 

LOPR  “  0  makes  R  ■  constant 

2 

LOPR  «  pos .  makes  R  ■  d/V  in  accordance  with  the 
Poi8euille-Hagen  law. 

PHEART  -  calculates  the  forcing  function  of  the  heart. 

PCHART  -  calculates  pressures  and  valve  dispositions  within  the 
heart . 

SKON  -  calculates  compliance  pressures.  Options  available  are 
L0PP1  -  0  means  compliance  constant.  L0PP1  +  0  makes  compliance 
a  function  of  volume.  (See  Appendix  1) 

HAMING  and  RUNGE  -  calculates  advancement  of  all  independent 


variables  in  time 
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Calling  structure: 


PHEART  RESIST 


BMAIN 

SKON  BDRIVS 

SKON  PCHART 


HAMING 


RUNGE 


BCC  -  program  of  gas  exchange  with  constant  blood  flow 
BCC  -  (same  as  BMAIN) 

Subroutines  of  BCC 

BCDR1V  -  calculates  derivatives  of  tissue  concentrations. 
Options  are: 

LOPGS  ■  1  execute  for  only  the  gases  0^  and  C0^ 

LOPGS  *  2  execute  for  only  the  gases  0^.  C0^»  and  X 

LOPGS  *  3  execute  for  only  the  gases  X 

LOPGS  ■  4  execute  for  only  the  gases  0^*  CO^,  and  CO 

LOPGS  -  5  execute  for  all  the  gases  0^,  C02,  CO,  and  X  where 

X  is  one  or  more  gases  that  do  not  bind  to  hemoglobin 


02DRIV  -  calculates  the 


dt  *  dt 


C02  „  dPC0  f  „ 

~rr —  ,  and  for  all 

dt  dt 


chambers.  (— rr—  is  not  calculated  if  LOPGS. LT. 4). 
dt 


XDRIV  -  calculates  the  -rr  for  all  chambers  (is  not  called  unless 

at 

LOPGS  -  2,  3,  or  5) 


dso2  dso2  dcco2  dcco2 

SATUR  -  calculates  Sn  ,  Crn  ,  -rs —  ,  -r= —  ,  -r= -  ,  and  -r= - 

°2  C02  dP02  dPCC>2  dP02  dPC02 


GFLTST  -  calculates  Kn  and  K__  for  all  chambers. 

°2  LU2 

FLTEST  -  determines  the  direction  of  flow  and  (Cv)in  and  (Cv)out 


for  each  chamber 


BCC 


Calling 


BCMAIN 


MAIN 


I 


Structure 


BCDRIV  Subroutines 

of 

BMAIN 

02DRIV  XDRIV 
SATOR  GFLTST  GFLTST 

FLTEST  FLTEST 

program  of  gas  exchange  with  pumping  heart  (circulatory 
model  with  gas  transport) 

BCMN  -  same  as  BMAIN. 

Subroutines  of  BCMAIN  consist  of  all  those  in  BMAIN  plus  all 
those  in  BCC. 

complete  cardio-pulmonary  model  to  date. 

MAIN  -  controls  the  execution  of  the  blood  program  (BCMAIN) 
and  lung  program  (LMAIN)  so  that  they  may  be  executed  with 
different  time  steps  and  durations. 

Subroutines  in  MAIN. 

BCMAIN  -  circulatory  model  with  gas  transport  with  various 
entry  points  added. 

LMAIN  -  complete  lung  program  with  various  entry  points 
added. 

BLCONV  -  calculates  the  rate  of  removal  of  each  gas  from 
each  alveoli  chamber. 

LBCONV  -  calculates  external  pressures  on  blood  chambers 
using  alveoli  pressures  and  pleural  pressures  from  the  lung 


mu 


■N.  \  V  \ 
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program  and  converts  each  gas  concentration  to  each  gas 
partial  pressure  in  each  alveoli. 


Given 


PART  II 

DEVELOPMENT  OF  A  NONINVASIVE  METHOD  FOR  DETERMINING 
LONG  VENTILATION  DISTRIBUTIONS  WITH  THE  AID  OF  A 
NUMERICAL  SIMULATOR  OF  THE  HUMAN  PULMONARY  SYSTEM 

I.  INTRODUCTION 

There  are  numerous  techniques  for  diagnosing  and  analyzing  lung 
diseases  which  cause  uneven  ventilation.  These  can  involve  a  wide 
variety  of  medical  diagnostic  techniques  ranging  from  the  common  chest 
x-ray  to  the  more  elaborate  technique  of  fiberoptic  bronchoscopy. 
This  report  presents  a  new  method  for  diagnosing  and  analyzing  lung 


diseases  which  cause  uneven  ventilation.  This  new  method  is  con¬ 
sidered  advantageous  relative  to  most  other  methods  for  three  main 
reasons:  it  is  noninvasive,  it  is  relatively  inexpensive  to  use,  and 
it  can  be  performed  in  a  relatively  short  period  of  time.  Before 
describing  this  new  method  some  of  the  other  methods  currently  used 
for  diagnosing  lung  disease  are  reviewed  below  for  later  comparison 
with  the  new  method  Introduced.  This  is  preceded  by  a  brief  explana¬ 
tion  of  pulmonary  ventilation  and  associated  diseases. 

The  most  important  function  of  the  lungs  is  to  ventilate  the 
alveoli,  the  functional  units  of  the  lung  parenchyma,  for  the  purpose 
of  gas  exchange  with  the  blood.  Thus,  tissues  of  the  body  are  sup¬ 
plied  with  oxygen  and  relieved  of  carbon  dioxide.  Lungs  also  elimi¬ 
nate  waste  in  vapor  or  gas  form  and  serve  as  a  blood  reservoir  for 
controlling  fluctuations  in  extrapulmonary  blood  volume  and  pressure. 
The  volume  ventilating  the  lungs  each  breath  is  approximately  equal  to 
the  volume,  known  as  the  tidal  volume  (TV),  entering  the  nose  and 
mouth  during  each  breath  with  some  correction  due  to  changes  in 
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temperature  and  water  vapor  pressure  as  air  enters  the  lungs.  If 
ideal  ventilation  were  to  occur  inspired  air  of  the  samd  composition 
would  be  distributed  uniformly  to  each  of  the  parenchymal  units,  each 
unit  receiving  an  equal  fraction  of  its  preinspiratory  volume.  In 
reality  ideal  ventilation  does  not  occur,  even  in  healthy  persons. 
However,  certain  environmental  factors  and  disease  states  can  cause 
extremes  in  the  imbalance  of  ventilation  of  different  parts  of  the 
lungs.  Asthma,  emphysema,  bronchitis,  and  tumors  cause  uneven  ven¬ 
tilation.  Changes  in  whole  body  acceleration,  such  as  that  experi¬ 
enced  by  astronauts  and  pilots  of  high  performance  aircraft,  frequent¬ 
ly  have  detrimental  effects  on  regional  pulmonary  ventilation  that  may 
contribute  to  impairment  of  abilities  to  function  physically  and 
mentally  (Burton,  et  al.,  1976).  In  all  cases  of  uneven  ventilation 
the  result  is  impaired  or  ceased  function  of  the  poorly  ventilated 
lung  regions. 

A  number  of  techniques  are  presently  used  for  analyzing  pulmonary 
ventilation.  In  clinics  and  pulmonary  function  laboratories  the 
techniques  primarily  used  include  the  inert  gas  washout  test,  the 
helium  distribution  test,  the  single  breath  oxygen  test,  and  bronchlo- 
splrometry  (Comroe,  et  al. ,  1962).  Each  of  these  techniques  may  be 
utilized  to  augment  an  examination  of  regional  lung  function.  How¬ 
ever,  of  these  methods  only  bronchlospirometry  provides  a  means  for 
analyzing  regional  ventilation,  and  this  is  only  a  gross  analysis.  It 
aids  only  in  distinguishing  the  ventilation  of  the  major  divisions  of 
the  lungs  by  catheterization  of  one  or  both  of  the  main  bronchi.  None 
of  these  methods  provides  a  means  of  analyzing  the  ventilation  charac¬ 
teristics  within  a  region  smaller  than  the  two  major  divisions. 


Some  of  the  techniques  used  in  research  provide  a  more  detailed 
assessment  of  regional  ventilation  of  the  lung.  The  volume  of  gas  in 
poorly  ventilated  or  unventilated  alveoli  can  be  calculated  by  using 
the  body  pie  thy  smo  graph  to  measure  the  total  volume  of  gas  in  the 
thorax  and  subtracting  from  this  the  volume  measured  by  a  technique 
that  washes  out  only  alveolar  gas  in  communication  with  airways 
(Comroe:  Physiology  of  Respiration,  1974).  A  commonly  employed 
method  of  determining  regional  ventilation  is  performed  by  measuring 
regional  clearances  of  inhaled  Xenon- 133,  or  other  radioactive  gas, 
using  externally  positioned  scintillation  detectors  (Schmid,  et  al., 
Jap,  49:279-286,  1980;  Rehder,  et  al.,  Jap,  42:391-402,  1977; 
Dolfuss,  et  al.,  Resplr.  Physiol.,  2:234,  1967).  In  a  similar  method 
Xenon- 133  dissolved  in  saline  is  injected  intravenously,  and  the 
washout  of  this  poorly  soluble  gas  is  measured  with  external  detectors 
as  it  diffuses  into  the  alveoli  (Dollery  and  Hugh-Jones,  Brit.  Med. 
Bull.  19:59,  1963).  Another  more  accurate  method  is  the  use  of 
fiberoptic  bronchoscopy  arranged  for  sampling  gasses  at  specified 
locations  within  the  bronchial  tree  (Denison,  et  al. ,  Resp.  Physiol. 
42:87-99,  1980). 

The  research  methods  discussed  above  are  used  to  determine 
whether  uneven  ventilation  is  occurring  and  approximately  to  what 
extent  it  is  occurring.  However,  of  these  methods  only  the  radio- 
graphic  techniques  and  the  recently  tested  fiberoptic  bronchoscopy 
provide  information  regarding  the  locations  of  the  poorly  ventilated 
regions.  Experiments  with  various  research  laboratory  techniques 
occasionally  provide  new  approaches  for  estimating  regional  lung 
ventilation,  but  few  of  the  research  techniques  currently  used  are 


wAm 


Parallel  Ventilation 

(D.  M.  Gomez,  RN.A.S.  49:312,  1963) 


Figure  1  Schematic  of  the  Gomez  Ventilation  Model 
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rapid  methods  employing  safe,  noninvasive  procedures  appropriate  for 
clinical  application. 

The  initial  objective  of  this  work  is  to  introduce  a  noninvasive 
technique  for  the  indirect  determination  of  the  distribution  of 
ventilation  within  the  lungs  from  gas  concentrations  in  expired  air 
measured  at  the  mouth.  It  is  believed  that  such  a  technique  can  be 
further  developed  as  a  noninvasive  procedure  for  determining  regional 
ventilation. 

The  idea  that  changes  in  the  distribution  of  ventilation  within 
the  lungs  can  be  determined  from  the  record  of  changes  in  concen¬ 
trations  of  test  gas  in  expired  air  measured  at  the  mouth  is  not  a  new 
concept.  Domingo  M.  Gomez  proposed  such  a  technique  in  previous 
reports  (Gomez,  1962,  1963),  but  because  of  limitations  of  his  method 
it  has  not  been  widely  used.  Gomez  proposed  treating  the  lungs  as  a 
continuous  collection  of  volume  units  with  parallel  ventilation,  with 
each  unit  having  some  value  for  ventilation,  Z,  defined  as  the  ratio 
of  tidal  volume  to  functional  residual  capacity.  This  model  is  shown 
schematically  in  Figure  1.  Using  this  assumption  Gomez  showed  that, 
for  washout  of  a  non-diffusing  gas,  the  graph  of  expired  concen¬ 
tration,  C (N) ,  versus  number  of  breaths,  N,  was  uniquely  determined  by 

* 

the  distribution,  F(Z),  of  ventilation  in  the  lungs.  In  other  words 
each  volume  element  of  the  lungs  has  a  value  for  the  ratio  Z,  and 
F(Z)DZ  is  the  fraction  of  all  chambers  having  values  for  this  ratio 
between  Z  and  Z+DZ.  Gomez  then  showed  that  one  could  determine  the 
graph  of  F(Z)  versus  Z  for  a  subject  from  the  graph  of  C(N)  versus  N 
obtained  in  an  inert  gas  washout  test  as  shown  in  Figure  2, 
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Figure  2  Calculatlonel  Scheme  in  the  Gomez  Model 
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The  Gomez  method  of  calculating  F(Z),  ventilation  distribution, 
from  C(N),  expired  gas  concentrations  measured  at  the  mouth,  has  been 
used  to  compute  the  distribution  of  ventilation  in  the  lung  even 
though  the  actual  lung  is  not  ventilated  exclusively  in  parallel,  as 
in  the  Gomez  model.  Thus  the  computed  distribution,  F(Z),  is  the 
equivalent  distribution  that  parallel  ventilated  lungs  would  need  in 
order  to  function  as  the  actual  system  with  respect  to  ventilation. 
Denison,  et  al.  (Resp.  Physiol.  42:87,  1980)  used  information  from 
lung  gas  samples  taken  from  inside  the  lung  via  bronchoscopy  to 
demonstrate  changes  in  regional  ventilation  resulting  from  changes  in 
body  position  as  indicated  in  Figure  3.  Hence,  changes  in  the  graph 
of  F(Z)  versus  Z  with  change  in  body  position  may  be  interpretable. 
In  fact,  it  may  be  possible  to  determine  not  only  the  fractional 
distribution  of  ventilation  but  also  to  determine  which  spatial 
regions  have  large  or  small  ventilation  solely  by  this  noninvasive 
method  simply  by  examining  changes  in  concentrations  of  expired  gasses 
resulting  from  changes  in  body  position. 

The  suggested  procedure  is  to  compute  the  ventilation  dis¬ 
tribution,  F(Z),  from  expired  gas  concentration,  C(N),  at  a  given  body 
position.  Then  change  the  body  position  and  repeat  the  run.  This  is 
repeated  for  a  number  of  body  positions  to  produce  a  family  of  ven¬ 
tilation  distribution  curves  computed  for  a  particular  subject.  It  is 
anticipated  that  changes  in  the  general  pattern  formed  by  this  family 
of  curves  can  be  correlated  with  changes  in  regional  ventilation  in 
the  lungs,  as  indicated  in  Figure  4. 

A  computerized  simulator  for  the  human  lungs  is  being  used  for 
this  research.  A  number  of  mathematical  simulators  have  been 
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Figure  3  Conceptualization  of  Denison  et  al.  Ventilation  Test 
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II.  Application  of  Technique 


Figure  4  Conceptualisation  of  Method  for  Determining  Local 
Ventilation  Dya functions 


developed  for  modeling  the  lungs  (Collins.  1967;  Grimby,  et  al., 
1968;  Saidel,  et  al..  1971;  Lambert  and  Wilson,  1973;  Fry  and 
Hyatt,  1960;  West,  1969;  Mead,  et  al.,  1970;  Saidel,  et  al.,  1972; 
and  Hilberman,  et  al.,  1973).  One  that  serves  as  a  more  complete 
numerical  simulator  of  the  human  pulmonary  system  has  been  recently 
developed  (Calvert,  1980).  A  version  of  this  computerized  model  has 
been  used  to  provide  inert  gas  washout  data  for  this  research.  The 
model  provides  a  convenient  means  for  performing  experimental  trials 
for  this  investigation,  because  known  variations  in  regional  ventila¬ 
tion  can  be  programmed  to  the  model  as  well  as  changes  in  body  posi¬ 
tion  and  whole  body  acceleration.  Clinical  tests  have  yet  to  be 
performed  to  provide  inert  gas  washout  data  from  subjects  with  known 
regional  lung  characteristics  to  verify  the  methods.  This  concept  is 
illustrated  in  Figure  5. 

The  research  plan  included  the  following: 

1)  The  algorithm  of  the  computerized  numerical  simulator  of  the 
lungs  was  amended  so  that  it  provided  data  simulating  a  multi¬ 
breath  inert  gas  washout  test. 

2)  A  companion  program  named  N2WASH  was  written  that  retains  values 
for  breath  order  N  and  mean  concentration  C(N)  of  inert  gas  in 
expired  air  of  each  breath.  This  program  also  uses  data  from  the 
lung  model  to  determine  the  tidal  volume,  TV,  and  functional 
residual  capacity,  FRC,  of  each  breath  during  the  simulated  inert 
gas  washout. 

3)  The  best  fit  for  a  curve  of  C(N)  versus  N  was  computed,  and  the 

data  from  this  curve  was  used  to  compute  the  distribution,  F(Z) , 

of  ventilation,  TV/FRC.  At  this  point  in  the  investigation  a 
flaw  was  discovered  in  the  numerical  methods  used  to  compute  the 
ventilation  distributions  using  parameters  computed  from  the 
numerical  approximation  of  the  nitrogen  washout  curve.  (An 

explanation  of  this  flaw  in  the  methods  is  provided  below.)  The 

body  position  simulated  in  the  lung  model  was  changed  and  the 
washout  test  repeated.  Computation  of  a  distribution  was  planned 
for  each  body  position,  and  thus  a  profile  of  distributions  would 
be  established  to  characterize  the  lungs  simulated. 
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4)  Step  3  vas  to  be  performed  with  simulated  normal  lungs  and  with 
regional  ventilation  dysfunctions  programmed  into  the  lung  model. 
The  completed  procedure  was  expected  to  provide  -  a  separate 
profile  of  distributions,  F(Z),  representing  each  different 
ventilation  dysfunction.  The  most  critical  part  of  this  step 
would  be  to  determine  a  relation  between  the  computed  dis¬ 
tributions  F(Z)  and  corresponding  regional  ventilation  patterns 
that  allow  accurate  prediction  of  these  patterns. 

It  was  discovered  that  the  flaw  in  the  numerical  methods  lay  in 
the  equation  used  to  approximate  the  nitrogen  washout  curve.  This 
equation,  a  sum  of  exponential  terms,  was  sensitive  to  error  in 
initial  approximations  of  the  curve  fitting  parameters.  In  order  to 
obtain  parameters  which  provided  a  reasonably  close  approximation  to 
the  expected  ventilation  distribution,  the  initial  estimates  of 
parameters  for  the  curve-fitting  equation  had  to  be  very  close  to  the 
final  values  computed  with  an  iterative  least  squares  procedure. 
Otherwise  the  computatibns  failed  to  close  in  on  acceptable  parameter 
values,  or  a  large  and  costly  number  of  iterations  was  required  to 
close  in  on  acceptable  values  for  the  parameters.  The  flawed  itera¬ 
tive  procedure  was  the  application  of  the  afore-mentioned  Lagrangian 
method  of  undetermined  multipliers  used  for  minimization  in  a  least 
squares  procedure.  During  each  iteration  new  values  were  selected  for 
the  curve  fitting  parameters  using  an  interval  halving  approach. 

Recently  an  alternative  method  for  approximating  lung  ventilation 
distribution  curves  has  been  developed.  This  alternative  procedure  is 
described  below  in  a  subsequent  portion  of  this  report.  At  present 
this  alternative  procedure  is  being  implemented  for  achieving  the 
ultimate  goal  of  this  investigation,  that  of  developing  a  noninvaslve 
procedure  for  determining  the  locations  of  localized  or  regional 
dysfunctions  within  the  lungs. 
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Development  of  a  nonlnvaslve  method  of  characterizing  ventilation 
distributions  could  lead  to  a  useful  tool  for  clinical  applications. 
In  particular  this  method  may  be  useful  in  diagnosing  regional  ven¬ 
tilation  dysfunctions  such  as  occur  with  emphysema,  fibrosis  and 
tumors . 

II.  EQUATIONS  AND  COMPUTATIONS  METHODS 

A.  Mathematical  Methods 

The  initial  goal  of  this  research,  as  stated  in  the  introduction, 
is  to  develop  a  numerical  method  for  determining  the  ventilation 
distribution  in  the  human  lungs.  The  suggested  procedure  requires  the 
measurement  of  end-expiratory  nitrogen  concentrations  taken  at  the 
mouth  during  a  nitrogen  washout  procedure.  The  next  step  in  determin¬ 
ing  the  ventilation  distribution  is  to  fit  the  nitrogen  washout  data 
with  an  equation  for  the  nitrogen  washout  curve.  This  curve-fitting 
is  performed  by  a  least  squares  procedure.  However,  this  procedure  is 
complicated  by  two  constraints  imposed  by  the  Gomez  equations  (Gomez, 
1963)  used  in  this  investigation  for  determining  the  ventilation 
distribution.  The  curve  fitting  equation  and  the  two  constraint 
equations  are  given  below. 

The  curve  to  fit  is: 

C(N)  -  VL/VT*[SUM  OVER  J-1,NPJ  A(J) ((1/((B(J,L)  -  2+N)**(M(J,L)+l))) 

-  <1/«B(J,L)  -  1+N)** (M(J ,L)+1))))  ,  N  -  1,2 . NBO.  [1] 

The  two  constraints  are: 


1  -  [SUM  OVER  J-l,  NP]  A(J)(l/((B(J,L)  -  1)**(M(J,L)+1)))  [2] 

and 

VT/VL  -  [SUM  OVER  J-l,  NP]  A(J) ((1/((B(J,L)  -  2)**(M(J,L)+1))) 

-  (1/((B(J.L)  -  1)**(M(J,L)+1))))  [3] 

The  two  constraint  equations  require  the  use  of  Lagrange's  method 
of  undetermined  multipliers  along  with  the  least  squares  procedure  to 
select  parameters  for  the  best  fitting  washout  curve.  The  parameters 
selected  in  this  fashion  are  given  values  which  are  useful  also  in 
constructing  the  ventilation  distribution  curve. 

The  least  squares  method,  applied  in  conjunction  with  Lagrange's 
method  for  unknown  multipliers  is  used  to  find  the  values  for  the 
parameters  A(J)  which  provide  the  best  fitting  approximation  for  the 
washout  curve.  This  least  squares  approach  generates  a  number  of 
simultaneous  equations  which  are  solved  later  in  a  subsequent  elimina¬ 
tion  routine.  The  coefficient  and  solution  matrices  used  in  the 
elimination  are  established  in  the  program  CRVFIT  (see  section  below 
on  program  CRVFIT)  as  follows: 
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where  NP  is  the  number  of  independent  variables  A(J) ,  and  NP+2  is  the 
total  number  of  variables  A(J)  Including  the  two  variables  (two 
Lagrangian  undetermined  multipliers)  determined  by  the  two  constraint 
equations.  The  D(I,J)  are  constants  derived  in  the  least  squares 
procedure.  The  first  two  rows  contain  coefficients  corresponding  to 
constraint  equations  and  should  be  placed  at  the  top  of  the  matrix, 
because  the  last  two  coefficients  in  these  rows  are  zero.  The 
elimination  routines  require  that  no  zeroes  lie  on  the  diagonal.  The 
constants  P(J)  represent  the  data  vector. 

Initial  values  for  parameters  M(J,L)  and  B(J,L)  must  be  input  by 
the  user.  The  minimum  dimension  required  for  subscript  J  in  both 
M(J,L)  and  B(J,L)  is  two  less  than  the  number  of  unknown  parameters 
A(J)  used  in  the  curve  fitting  process. 

The  equation  for  the  curve  to  fit,  shown  above,  is  derived  from 
mass  balance  equations  for  a  nitrogen  washout  procedure.  In  the 
process  of  constraining  the  mass  balance  equations  to  accommodate  the 
requirements  described  by  D.M.  Gomez  (Gomez,  1963)  for  determining  the 
ventilation  distribution,  the  following  conversion  was  made  for 
convenience  in  deriving  those  mass  balance  equations: 

U  -  LN(1+Z) 

Here  Z  represents  the  ratio  of  TV/FRC  for  a  ventilating  unit  in  the 
lungs.  Using  this  conversion  the  following  model  was  derived  for  the 
ventilation  distribution: 


F(Z)  -  F1(U) 
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F1(U)  -  [SOM  OVER  J]  (A(J)/M(J ,L! )*(U**M(J,L))*EXP(-B(J,L)U) 

Initial  values  for  the  constants  H  and  B  are  chosen  in  the 
following  nanner: 

The  constants  M(J,L)  and  B(J,L)  are  chosen  such  that  the  plot  of 
fractional  distribution  of  ventilation  F(Z)  versus  ventilation  Z  of 
healthy  normal  lungs  assumes  the  expected  narrow  unimodal  shape  (S.M. 
Lewis,  et  al.,  1978).  These  constants  can  be  selected  by  the  follow¬ 
ing  procedure:  a  point  is  selected  from  the  distribution  curve  where 
the  distribution  of  Z(Z  ■  TV/FRC)  is  expected  to  fall  to  a  negligible 
value.  For  example  it  may  be  that  Z  were  rarely  expected  to  be  10  or 
greater  under  experimental  conditions.  Therefore,  M  and  B  may  be 
chosen  such  that  the  value  for  F(Z)  at  the  point  where  Z  ■  10  is 
negligible  relative  to  the  maximum  value  for  F(Z).  The  maximum  value 
for  F(Z)  corresponds  to  the  Z  value  where  D/DZ  F(Z)  »  0  or  to  the  U 
value  where  D/DU  F1(U)  *  0.  Assuming  one  variable  term  A(J)  is  used, 

F(Z)  -  (A/M!)((LN(1+Z))**M)EXP(-B(LN(1+Z)))  or 
F1(U)  -  (A/M! ) * (U**M) *EXP (-BU) 

setting  D/DU  F1(U)  -  0 
yields  U  ■  M/B 

Therefore,  U  ■  M/B  at  the  point  of  maximum  distribution  F1(U).  Thus, 
for  negligible  F(Z)  at  Z  -  10: 


(F(Z)  at  Z  -  10] 


[F(Z)  at  max  (U  -  M/B)] 


or,  as  an  example,  let 


[F(Z)  at  Z  -  10]/[F(Z)  at  max  (U  -  M/B]  -  0.01 
(A/Ml ) (U**M) EXP (-BU) / (A/M! ) ( (M/B) **M) EXP (-M)  -  0.01 

this  yields 

EXP (M-BU) * ( (BU/M) **M)  -  0.01  . 

For  Z  -  10,  U  -  LN(ll) .  Thus  the  appropriate  values  for  M  (values  for 
M  which  are  expected  to  provide  the  narrow  unlmodal  shape  previously 
mentioned)  may  be  chosen  arbitrarily  and  the  corresponding  values  for 
B  found  by  trial  and  error. 

It  was  noted  that  values  chosen  for  the  constants  M(J)  and  B(J) 
should  be  selected  such  that  the  different  ratios  M(J)/B(J),  the 
ratios  representing  the  values  for  Z  where  F(Z)  is  at  a  maximum,  are 
different  enough  to  allow  for  the  sum  of  widely  spaced  curves  summed 
to  fit  the  data. 

B.  Mainframe  System  and  Language  Used 

The  computer  language  used  for  this  project  was  the  most  current 
version  of  Fortran  available  at  the  University  of  Texas,  Fortran  77. 
The  compiler  used  was  the  FTN5  compiler.  This  was  the  one  compiler 
available  at  the  university  providing  compiled  versions  of  Fortran 
programs  which  allowed  the  most  rapid  execution  on  the  mainframe 
computer  system.  The  mainframe  system  maintained  by  the  Computation 
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Center  at  the  university  was  a  Control  Data  Corporation  6600  Digital 
Computer. 


C.  The  Lung  Model 

***************  Program  LUNG  *************** 

Program  LUNG  is  a  numerical  simulator  of  the  human  pulmonary 


system  written  in  Fortran.  This  simulator  was  developed  by  Ray  E. 
Calvert  (1977,  1980)  in  his  research  efforts  for  both  the  Masters  and 
Doctoral  degrees  in  biophysics.  This  research  was  performed  under  the 
supervision  of  R.  Eugene  Collins,  Professor  of  Biophysics  at  the 
University  of  Houston.  the  model  provides  a  dynamic  mathematical 
representation  of  the  pulmonary  structures  and  their  functions  during 
respiration,  including  the  gross  pulmonary  structural  mechanics, 
mechanics  of  airflow  in  the  lung  airways,  and  ventilation  of  major 
lung  regions,  the  model  shown  schematically  in  Figure  6  is  composed 
of  six  alveolar  chambers  representing  the  alveoli  or  acini  ventilated 
through  five  airway  chambers  representing  the  bronchial  tree.  These 
airways  Interconnect  the  alveolar  chambers  in  an  arrangement  that 
provides  both  parallel  and  series  ventilation  of  the  alveolar  cham¬ 
bers.  In  addition  there  is  a  twelfth  chamber  representing  a  spiro¬ 
meter  which  is  used  as  the  sole  source  of  gasses  inspired.  the 
diaphragm  and  abdominal  and  thoracic  respiratory  muscles  are  modeled 
as  a  system  of  pistons  which  provide  the  driving  forces  of  inspiration 
and  expiration.  The  model  includes  the  fluid  pressure  within  the 
chest  and  abdominal  cavities  resulting  from  gravitational  forces  on 
the  thorax,  diaphragm  and  abdomen.  The  model  also  simulates  the 
hydrostatic  pressure  gradient  within  the  pleural  fluid  and  the  effects 
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Figure  6  Schematic  of  Ventilation  Compartments  in  Pulmonary  Simulator 
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it  has  upon  the  regional  lung  volumes.  Variations  in  the  magnitude  of 
the  gravitational  field  and  in  the  orientation  of  the  lungs  within  the 
field  are  simulated  by  the  model.  A  more  detailed  description  of  this 
model  is  given  in  previous  AFOSR  reports  on  this  project. 

The  main  program  coordinates  19  subroutines.  All  are  controlled 
by  22  input  options  which  allow  the  user  to  orchestrate  execution, 
input,  and  output  such  that  innumerable  combinations  of  lung  type, 
inspired  gasses,  and  breathing  dynamics  can  be  modeled. 

D.  Companion  Programs  for  Data  Processing 

***************  Program  LUNGID  *************** 

Program  LUNGID  is  another  program  written  by  Calvert  (1980)  to 
organize  data  output  by  program  LUNG.  This  program  reads  a  file  on 
which  data  from  program  LUNG  have  been  stored  and  organizes  and 
outputs  those  data  in  readable  form.  It  provides  large  volumes  of 
detailed  data  pertaining  to  a  great  number  of  the  mechanical  and 
quantitative  changes  in  the  lungs  that  occur  during  breathing.  This 
is  helpful  when  diagnosing  execution  errors. 

***************  Program  N2WASH  *************** 

This  program  reads  data  from  a  file  output  by  program  LUNG  and 
selects  and  prints  only  output  pertinent  to  nitrogen  washout  studies. 
Program  N2WASH  also  outputs  data  to  another  file  to  be  used  later  for 
the  purpose  of  curve-fitting  data  to  an  equation  for  the  nitrogen 
washout  curve  and  eventually  for  the  purpose  of  computing  ventilation 


distributions. 
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E.  The  Curve-Fitting  Programs 

***************  Program  CRVFIT  *************** 

This  program  selects  data  from  a  nitrogen  washout  procedure  and 
computes  a  best  fit  curve  for  the  nitrogen  washout  data.  The  equation 
used  to  construct  this  curve  is  an  exponential  equation  with  parame¬ 
ters  selected  to  conform  to  constraints  imposed  by  the  Gomez  equations 
(Gomez,  1963)  for  determining  the  ventllat.  1  distribution. 
Lagrange's  method  of  undetermined  multipliers  is  used  here  in  conjunc¬ 
tion  with  a  least  squares  approach  for  fitting  a  curve  with  two 
constraints.  Thus  parameters  selected  for  the  best  fitting  washout 
curve  possess  values  also  suitable  for  later  use  in  the  construction 
of  a  ventilation  distribution  curve.  As  explained  below  in  the 
section  on  program  VDPLOT  ventilation  distribution  curves  are  subse¬ 
quently  used  for  characterizing  the  lungs  tested.  The  washout  curve 
is  C(N)  versus  N  where  C(N)  represents  the  average  concentration  of 
nitrogen  in  the  lungs  at  the  end  of  breath  number  N.  The  equations 
for  the  curve  to  be  fitted  to  the  washout  data  are  repeated  here  for 
convenience : 


The  curve  to  fit  is: 


C(N)  -  VL/VT*[SUM  OVER  J-l.NP]  A(J)((1/((B(J,L)  -  2+N)**(M(J,L)+l))) 
-  <1/«B(J,L)  -  1+N)**(M(J,L)+1))))  ,  N  -  1,2 . NBO.  [1] 


The  two  constraints  are: 


1  -  [SOM  OVER  J-l,  NP]  A(J) ( 1/ ( (B(J,L)  -  1)**(M(J,L)+1)))  [2] 
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and 

VT/VL  -  [SUM  OVER  J-l,  NP]  A(J) ((1/((B(J,L)  -  2)**(M(J,L)+1))) 

-  (1/((B(J,L)  -  1)** (M(J ,L)+1))))  [3] 

An  iterative  procedure  has  been  developed  for  finding  the  optimal 
values  for  the  parameters  M(J,L),  B(JfL),  and  A(J)  for  the  nitrogen 
washout  curve  approximation.  The  following  procedure  is  used: 
initial  values  for  parameters  M  and  B  are  selected.  Using  these 
values  the  values  for  the  corresponding  A(J)  are  computed.  The 
corresponding  washout  curve  CEQ(L,N)  is  computed  and  compared  with  the 
data  values  C(N)  for  the  washout.  The  average  of  the  absolute  values 
of  the  deviations  of  the  computed  values,  CEQ(L,N),  from  the  observed 
values,  C(N) ,  is  computed  each  iteration,  this  average  deviation  is 
labelled  SD(J,L).  One  of  the  values  for  M  or  B  is  temporarily  incre¬ 
mented,  and  a  new  washout  curve  and  SD(J,L)  are  computed.  If  the 
SD(J.L)  Increases,  indicating  an  increase  in  error,  a  change  in  the 
most  recently  altered  parameter,  M  or  B,  is  selected  to  be  -%  times 
the  previous  change,  l.e., 

DM(J.L)  -  -J*  *  (DM(J,L  -  1) 


or 


DB(J,L)  -  -h  *  (DB(J,L  -  1) 


However,  if  the  SD(J,L)  decreases  initially  a  change  in  the  parameter 
M  or  B  is  selected  to  be  twice  the  previous  change: 

DM(J,L)  -  2.0  *  DM(J,L  -  1) 
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or 

DB(J,L)  -  2.0  *  DB(J,L  -  1) 

Use  of  this  interval-doubling  process  is  continued  with  each  parameter 
for  which  incrementation  has  resulted  in  a  decrease  in  average  de¬ 
viation,  SD(J,L).  After  an  incrementation  results  in  an  increase  in 
average  deviation,  from  that  time  onward  an  interval-halving  routine 
is  employed  until  the  desired  accuracy  is  achieved  or  until  the  time 
limit  provided  by  the  user  is  exceeded.  If  the  SD(J,L)  decreases 
while  using  this  routine  without  interval-doubling,  the  subsequent 
change  in  the  parameter  M  or  B  is  selected  to  be  Identical  to  the 
previous  change: 

DM(J.L)  -  DM(J,L  -  1) 
or 

DB(J,L)  -  DB(J,L  -  1) 

This  process  of  choosing  parameters  such  that  the  average  deviation, 
SD(J,L),  decreases  is  continued  until  the  average  deviation  is  less 
than  or  equal  to  an  arbitrarily  selected  maximum  called  SDMAX.  When 
the  average  deviation  becomes  smaller  than  this  arbitrarily  chosen 
value  the  approximated  washout  curve  is  considered  to  be  an  adequately 
accurate  fit  and  may  be  used  in  obtaining  information  about  ventila¬ 
tion  distribution.  This  information  is  subsequently  used  to  extract 
Information  concerning  regional  function  in  the  lungs. 
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***************  Program  FAMZ  *************** 

Program  FAMZ  has  two  different  uses.  One  is  to  plot  data  on  a 
Diablo  Teletype  terminal.  Thus  clinical  washout  data  or  modeled 
washout  data  may  be  compared  graphically  with  computed  washout  data 
provided  from  program  CRVFIT.  These  data  are  compared  by  plotting  all 
curves  on  the  same  graph.  This  is  convenient  when  the  user  wishes  to 
obtain  a  graphical  comparison  for  visualizing  the  closeness  of  fit  of 
an  approximated  curve  to  the  curve  being  fitted. 

The  second  use  of  Program  FAMZ  is  to  compute  and  plot  washout 
concentrations  for  up  to  five  different  washout  curves  on  a  single 
graph.  These  curves  are  computed  within  program  FAMZ  from  the  washout 
curve  equation  with  parameters  input  by  the  user.  This  use  of  program 
FAMZ  provides  a  tool  by  which  one  may  familiarize  oneself  with  the 
graphing  characteristics  of  the  equations  used  to  approximate  the 
changes  in  expired  nitrogen  concentration  observed  in  a  nitrogen 
washout  procedure.  These  equations  are  provided  in  the  section  above 
on  program  CRVFIT. 

The  approach  to  familiarization  is  to  execute  this  program  a 
number  of  times  changing  the  value  of  one  of  the  parameters  M(J)  and 
B(J)  before  each  execution  while  keeping  other  parameters  constant. 
According  to  the  curve  fitting  equations  the  A(J)  depend  upon  the  M(J) 
and  B(J)  chosen.  During  each  execution  a  fitted  washout  curve  of 
CEQ(L,N)  versus  N  is  plotted  to  provide  a  visual  means  by  which  to 
judge  the  effects  of  the  changes  employed,  this  is  done  with  the 
Hyplot  Plotting  Routine  provided  as  a  library  program  by  the  Universi¬ 
ty  of  Texas  for  creating  graphs  with  a  Xerox  Diablo  Computer  terminal. 
These  graphs  using  arbitrarily  chosen  parameters  may  be  compared  to 
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graphs  of  clinical  washout  data  from  normal  healthy  individuals  or 
graphs  of  data  from  the  lung  program  to  determine  what  values  for  the 
curve-fitting  parameters  may  be  considered  appropriate  for  modeling 
the  washout  of  normal  lungs. 

This  is  not  an  accurate  test  of  the  effects  of  the  parameters 
M(J)  and  B(J)  computed  by  program  CRVFIT,  because  the  A(J)  are  not 
accurately  computed.  The  A(J)  are  kept  constant  for  each  set  of 
parameters.  This  allows  investigation  of  the  effects  of  the  varying 
parameters  M(J)  and  B(J)  on  the  plots  without  interference  from 
changing  A(J).  To  visualise  the  washout  curve  resulting  from  more 
accurately  computed  parameters  from  program  CRVFIT,  all  three  sets  of 
parameters,  M(J),  B<J),  and  A(J)  are  input  by  the  user.  To  do  this  a 
read  statement  is  used  in  program  FAMZ  to  input  the  A(J),  and  the  code 
used  to  compute  the  A(J)  is  by-passed. 

F.  Computing  the  Ventilating  Distribution 
***************  Program  VENDIS  *************** 

This  program  reads  the  washout  curve-fitting  parameters,  M(J) , 
B(J),  and  A(J),  from  an  input  file  created  in  program  CRVFIT  and  uses 
the  ventilation  distribution  equation  to  compute  the  ventilation 
distribution  for  the  particular  lungs  which  were  used  to  obtain  the 
curve-fitting  parameters.  The  equation  for  ventilation  distribution 
is 


F(Z)  -  [SOM  OVER  J-l,  NP]  A(J)/M(J)I  *  [(LN(1+Z))**M(J)]/(1+Z)**B(J) 


Finally  the  arbitrarily  chosen  data  points  for  ventilation,  Z, 
and  the  corresponding  ventilation  distributions,  F(Z),  computed  in 
program  VENDIS  are  stored  on  an  output  file  for  later  use  in  the 
plotting  routine,  VDPLOT. 

This  program  uses  the  1MSL  routine  gammas,  which  computes  the 
gamna  function  for  given  input  values.  The  gamma  function  is  used  to 
compute  the  factorial  function  of  the  parameters  M(J) . 

***************  Program  VDPLOT  *************** 

This  program  reads  and  plots  ventilation  distribution  data  using 
the  Hyplot  routine  for  plotting  graphs  with  a  Diablo  terminal.  There 
are  as  many  as  24  sets  of  washout  data  for  each  individual  set  of 
lungs  examined.  These  correspond  to  the  24  different  seatback  angles 
used  for  the  washout  test.  For  a  simultaneous  plot  of  all  ventilation 
distributions  for  a  particular  set  of  lungs  this  program  requires  that 
data  files  for  all  seat  tilts  for  a  particular  set  of  lungs,  up  to  24 
files,  be  copied  to  a  single  file.  The  value  for  NPLOTS,  the  number 
of  graphs  to  be  plotted,  must  also  be  entered  on  the  input  file. 

The  programming  Included  in  this  routine  alters  input  data  such 
that  each  distribution  curve  is  slightly  offset  from  the  previous  one. 
This  provides  a  three-dimensional  effect  or  a  perspective  view  of  the 
curves.  All  the  curves  thus  plotted  together  can  be  visualized  as  the 
framework  of  a  surface.  The  appearance  of  this  framework  may  be  used 
to  suggest  the  physical  characteristics  of  the  lungs  it  represents, 
especially  uneven  ventilation  within  different  regions  of  the  lungs. 
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G.  Alternative  Method  for  Computing  Ventilation  Distribution 
A  fifty-compartment,  numerical  model  employed  by  Vagner,  et  al. 
(1974)  to  compute  distributions  of  ventilation-perfusion  ratios 
provides  an  algorithm  adaptable  for  computing  ventilation  dis¬ 
tributions.  This  method  used  by  Vagner  and  co-workers  computes  the 
distribution  of  ventilation-perfusion  ratios  from  a  data  curve  relat¬ 
ing  the  solubilities  and  the  changes  in  arterial  and  expired  air 
concentrations  of  several  inert  gasses  artificially  introduced  into 
the  bloodstream.  the  adapted  algorithm  computes  ventilation  dis¬ 
tributions  using  data  from  an  inert  gas  washout  curve.  Lewis  and  his 
colleagues  (1978)  have  used  a  similar  50-compartment  model  to  compute 
ventilation  distributions  from  inert  gas  washout  data  and  have  used 
the  results  to  identify  ventilation  distributions  which  correlate  with 
certain  lung  diseases.  A  similar  compartmental  lung  model  is  used  in 
this  investigation  to  compute  ventilation  distributions.  Instead  of 
using  clinical  data  to  compute  the  ventilation  distribution  curves  the 
lung  model  developed  by  Calvert  (1980)  is  used  to  provide  the  inert 
gas  washout  data.  The  50-compartment  model  is  not  yet  in  use  for 
computing  the  ventilation  distribution  due  to  some  difficulties  in 
adapting  the  least  squares  approach  used  to  analyze  the  inert  gas 


washout  data 
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